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ABSTRACT 

Analytical and numerical solutions are obtained for the equation of radiative transfer in ultra-relativistic 
opaque jets. The solution describes the initial trapping of radiation, its adiabatic cooling, and the transition 
to transparency. Two opposite regimes are examined: (1) Matter-dominated outflow. Surprisingly, radiation 
develops enormous anisotropy in the fluid frame before decoupling from the fluid. The radiation is strongly 
polarized. (2) Radiation-dominated outflow. The transfer occurs as if radiation propagated in vacuum, preserv- 
ing the angular distribution and the blackbody shape of the spectrum. The escaping radiation has a blackbody 
spectrum if (and only if) the outflow energy is dominated by radiation up to the photospheric radius. 
Subject headings: radiative transfer — relativistic processes — scattering — gamma-ray burst: general 



1. INTRODUCTION 

Powerful jets from compact objects can have significant op- 
tical depth to scattering. The foremost example is gamma- 
ray bursts (GRBs). They are emitted by hot ultra-relativistic 
outflows that remain opaque until they travel a large distance 
from the central engine. Where the jet becomes transparent, 
the trapped radiation is released and contributes to the GRB. 
Its spectrum is expected to be nonthermal because of dissipa- 
tive processes in the subphotospheric region. 

This "photospheric emission" is likely the main component 
of observed GRBs. Recent work provides significant support 
for this picture. Three heating mechanisms have been pro- 
posed to shape the photospheric spectrum: (1) internal shocks, 
(2) dissipation of magnetic energy and excited plasma waves 
(Thompson 1994; Spruit, Daigne, & Drenkhahn2001; Ioka et 
al. 2007), and (3) collisional dissipation (Beloborodov 2010; 
hereafter B10). The latter mechanism is straightforward to 
model from first principles and turns out to reproduce the 
canonical GRB spectrum with no fine-tuning of parameters 
(B10; Vurm, Beloborodov, & Poutanen 201 1). 

Modeling emission from opaque jets requires simulations 
of radiative transfer. Two methods have been developed for 
such simulations to date. First, solving the kinetic equations 
for the electrons and photons that interact via Compton scat- 
tering (Pe'er & Waxman 2005; Vurm et al. 2011). Second, 
tracking a large number of photons that propagate and (ran- 
domly) scatter in the jet (Giannios 2006; B10). Pe'er (2008) 
used an analytic approach and Monte-Carlo simulations to 
study individual short pulses of photospheric emission. 

None of these works attempted to use the standard trans- 
fer equation. This approach is developed in the present paper. 
The extension of radiative transfer theory to relativistic out- 
flows is straightforward (e.g. Castor 1972; Mihalas 1980). In 
Section 2, we write down the transfer equation that is well- 
behaved (and simplifies) in the ultra-relativistic regime. Then 
we solve this equation for isotropic (Section 3) and exact 
(Section 4) models of electron scattering. In parallel, we ap- 
ply the independent Monte-Carlo technique and compare the 
results. 

In Section 5 we separately consider the case where radia- 
tion dominates the outflow energy up to the photosphere. This 

1 Also at Astro-Space Center of Lebedev Physical Institute, Profsojuz- 
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regime is of interest for GRBs with extremely low baryon 
loading, as described by Paczyhski (1986) and Goodman 
(1986). In contrast to their expectations, the transfer near the 
photosphere is not complicated and has a simple analytical 
solution. 

2. TRANSFER EQUATION IN THE ULTRA-RELATIVISTIC REGIME 
2.1. Formulation of the problem 

We are interested in outflows with Lorentz factors Y ^> 1 
and velocities j3 = v/c 1. Below we consider radiative 
transfer in outflows that are steady and spherically symmetric. 
These assumptions are not restrictive in the ultra-relativistic 
regime, as discussed in Section 6, — even strongly variable 
and beamed jets may be described by this model. 

The transfer equation is well-behaved in the limit /3 — > 1 
when it is formulated for radiation intensity in the fluid frame 
(see Appendix A). This frame is comoving with the outflow 
at any radius r. Hereafter intensity is denoted by J„(r, (i, v) 
where v is the photon frequency, [i = cos 6, and 9 is the pho- 
ton angle with respect to the radial direction. The quantities 
v, /i, and /„ are mea sured in the fluid frame. When (3 — > 1, 
the transfer equation (1A1 lb simplifies to 
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The quantity S u appearing in equation (03 is the source func- 
tion in the fluid frame. It is determined by how radiation in- 
teracts with the outflow. For example, the simplest model of 
isotropic and coherent scattering gives S u — (1/2) J l v d[i 
(e.g. Chandrasekhar 1960). 

The quantity r„ in equation ([TJ approximately represents 
the outflow optical depth (cf. Appendix B). It is defined as 
T u = K u r/Y where k v is the opacity in the fluid frame. In 
GRB jets, electron/positron scattering strongly dominates the 
opacity around the spectral peak of the burst. The scattering 
opacity is k u = an where a is the scattering cross section and 
n is the proper density of the flow. It may be expressed in 
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terms of the rate of e outflow through the sphere Airr 2 in the 
lab frame, N e = Airr 2 n T f3c. Then t v is given by 



aN„ 
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When the outflow carries no positrons, N e remains constant 
with radius. It also remains constant if positron creation is 
balanced by annihilation (this situation takes place in colli- 
sionally heated jets, see BIO). 

Besides the simpler form of the transfer equation, the ultra- 
relativistic regime implies a principal change in the formu- 
lation of the transfer problem. Radiation in the lab frame 
is strongly collimated and essentially all photons stream out- 
ward. Inward radial motion in the lab frame requires p, < —13 
in the fluid frame, which corresponds to a small solid angle 
Aft = 27r(l — At any given radius, only a small fraction 
£>(T~ 2 ) of all photons have p < As long as T > 1, 

one can exclude the solid angle p, < —(3 from the transfer 
problem, i.e. neglect its contribution 0(T~ 2 ) when calcu- 
lating the source function S v , and solve the problem in the 
domain — (3 < p < 1. Then the outer boundary condition is 
not needed, as information cannot propagate from larger r to 
small r. 

Thus, only an inner boundary condition should be speci- 
fied for the ultra-relativistic transfer problem. If it is given 
at a sufficiently small radius r; n , the radiation may be as- 
sumed isotropic in the fluid frame (as demonstrated by the 
solutions presented below, radiation maintains isotropy where 
t v > 100). 

Equation CO uses radius r as an independent variable. 
Alternatively, it can be written in terms of the comoving- 
observer time t(r) which is related to r by dt — dr//3cT. 
Then the transfer problem takes the form of an initial-value 
problem. One can think of ultra-relativistic transfer as the 
evolution of intensity I(p), as seen by the comoving observer. 
This view is valid as long as radiation is not allowed to stream 
backward in time t (i.e. backward in r). The model is exact 
in the limit T — > oo (then the entire region — 1 < p, < 1 
is included in the transfer domain without violating causal- 
ity). In practice, equation dTJ and the neglect of photons with 
p < — 3 gives an excellent approximation when V > 10. The 
approximation is perfect for GRB jets. 

2.2. Transfer of energy 

The flow of radiation energy is described by the frequency- 
integrated intensity, 

rOO r OO 

/(/x, r) = l v dv, S(p, r) = S v dv. (4) 
Jo Jo 

Integration of equation (fl~|i over v gives the equation for 

dl , 9N dl , (S-I) 

where r(r) = T~ 1 K(r)r is defined using the effective opacity, 

K = I~ j K u I w dv. 

2 This fraction equals A£7/47r (2r) — 2 for a very opaque flow, t v — ¥ 
oo, where radiation isotropy is maintained in the fluid frame. Deviations from 
isotropy that develop with decreasing t„ make this fraction even smaller, see 
Section 3. 



Multiplying equation (O by 1 + /i and integrating over /i, 
one gets 



d(I + h) 
d\nr 



= -A(I + h)+g(I + 2h+I 2 ) + T(S -I ), 

(6) 



where 7 m (to = 0, 1, 2) are the moments of intensity, 
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I m (r) = - j ^,r)v m dfx, 



(7) 



and So(r) is the zero moment of the source function. The 
quantities 47rJ m are the components of the stress-energy ten- 
sor of radiation, and equation © in essence expresses the first 
law of thermodynamics for radiation (e.g. Castor 1972). 

In the case of coherent scattering (e.g. Thomson scattering 
by a plasma with a small Kompaneets' y-parameter) So = Iq- 
Then the last term in equation (|6]l vanishes — there is no heat 
exchange between the plasma and radiation, and the evolution 
of radiation with r is adiabatic. Radiation does PdV work 
and gradually gives away energy to the bulk kinetic energy of 
the outflow as it propagates to larger radii. This process of 
adiabatic cooling is easily described at large optical depths, 
t> 1, where radiation is nearly isotropic in the fluid frame, 
I\ = and I2 = Iq/S- Then equation © gives 



din Jp _ ! 4 
d\nr 3 ^ 
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This equation reproduces the law of adiabatic cooling of ra- 
diation (adiabatic index 4/3) in expanding volume. For ex- 
ample, if r(r) = const (g = 1, see eq. |2| volume ex- 
pands as r~ 2 and radiation energy density in the fluid frame 
U = 47r/o / c decreases as r~ 8 / 3 . The energy flux through the 
sphere 47rr 2 measured in the lab frame decreases as r~ 2 / 3 . 
Note that equation ^ is valid only if radiation is isotropic in 
the fluid frame. Strong deviations from isotropy (which occur 
at optical depths r 10 as shown below) change the rate of 
adiabatic cooling. 

The fact that radiation does work on the outflow implies that 
the outflow accelerates and g in equation (Q]l is not an inde- 
pendent parameter of the transfer problem. It may be treated 
as a given fixed parameter only if the outflow inertia is large 
enough, so that it cannot be significantly accelerated by radi- 
ation. This condition reads pc 2 3> U where p is the rest-mass 
density of the outflow. This condition will be assumed in Sec- 
tions 3 and 4, where we solve the transfer equation with g w 1 
[i.e. r(r) ~ const]. In Section 5, we consider the opposite 
regime and discuss the coupled dynamics of the outflow and 
radiation. 

2.3. Transfer of photon number 

The photon number intensity is described by the following 
quantity, 



l(p,r) = 



■ dv, 



(9) 



where h is Planck constant. Integration of equation ([]]) over 
d In v gives the equation for X(/x, r), 



dl 



dl 



dlnr- (l-^^-Sd-wJJ+r^, (10) 

where S — §(S v /hv)dv and r(r) = T _1 K(r)r is defined 
using k = l~ l J n v Xdv. 
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Scattering conserves photon number, and the transfer equa- 
tion is expected to give the corresponding conservation law. 
Multiplying equation ( TTOb by 1 + \i, integrating over \x, and 
re-arranging terms, one gets 



d In i 



ln[r(J 1 +I )] = -2 ) 



(11) 



where T m (r) are the moments of I(/_t, r) and we used To = 
So which is true for any scattering process. The quantity Airli 
is the number flux of photons measured in the fluid frame, and 
47r2o /c is the number density of photons in the fluid frame. 
Lorentz transformation of the four-flux vector Att(Iq , X\ , 0, 0) 
gives the radial photon flux measured in the lab frame, T = 
AwT(Ii +j3Xo), Equation (fTTT i in essence states r 2 T — const 
(with /3 —> 1) and expresses conservation of photon number. 

3. ISOTROPIC-SCATTERING MODEL 

In this section, we solve the transfer problem assuming the 
simplest form of the interaction between radiation and the 
fluid: coherent isotropic scattering in the fluid frame. It gives 
a reasonable first approximation to Thomson scattering that is 
considered in Section 4. We consider here matter-dominated 
outflows — the outflow is assumed to be massive enough, so 
that it can coast with T(r) w const (Section 2.2), which cor- 
responds to g(r) m 1 (eq.[2j. 

Then the energy transfer equation (O reads, 



81 
8 In 7 
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Here we substituted the source function that describes 
isotropic scattering S(u,r) = Io(r), where Io is the zero- 
moment of intensity (eq. |7J- W e w iH assume a constant cross 
sectiorQ a(u) = const and N e (r) = const. Then equa- 
tion (f3]i gives 



t{t) = — , 
r 



R* = 
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Transfer of photon number is described by equation similar 
to equation (fl2l) where I is replaced by I and the numerical 
coefficient —4 in the second term on the right-hand side is 
replaced by —3 (cf. eq.fTOii. 

3.1. Integration of transfer equation 

Equation (fT2l gives the expression for 81/8 In r in terms 
of I. Direct integration in In r immediately yields the solu- 
tion for I(fi, r). Our numerical integration starts at r ln = 
3 x 10 -3 i?* and takes the isotropic I(u, r; n ) = const as the 
boundary condition. We use a uniform grid in 9 and In r of 
size 300 x 10 5 . With a simplest integrator — Runge-Kutta 
scheme of fourth order — the grid gives excellent accuracy 
of ~ 0.1% (we have checked this by varying the grid). Two 
more details of numerical integration are worth mentioning: 

(1) At one boundary of the computational domain \i — y — 1 
and the transfer equation gives (S — I) /I — > 0. This re- 
quires I = S = Jo at /j = —1. Note that the optical 

3 This is a good approximation for the bulk of GRB photons. The typical 
energy of observed photons is ~ 1 MeV. They are emitted in the rest frame 
of the jet with energy ~ r —1 MeV, much smaller than m e c 2 . Klein-Nishina 
corrections are small for such photons and the scattering cross section is ap- 
proximately independent of v. 




FIG. 1. — Angular distribution of radiation intensity in the fluid frame at 
three radii: r/R* m 0.03, 0.1 and 0.5, which correspond to optical depths 
t ss 30, 10, and 2, respectively. For a better comparison, we plot (r/R±) 2 I, 
where the factor (r/i?*) 2 compensates for the photon dilution due to ex- 
pansion. The overall normalization of the transfer solution is chosen so that 
I(fj,,r) = (r/i?*) -8 / 3 = t 8 / 3 at radii r < 0.01-R* where radiation is 
nearly isotropic (/ does not depend on fi) and follows the adiabatic cooling 
law / = T 8 / 3 . Open circles show the intensity at r = 2 that is obtained by 
the Monte-Carlo code (Section 3.3). 

depth Ar ray passed along the ray in one step A In r depends 
on fi: Ar ray (/i) — r(r) A In r/(l + /i). Numerical inte- 
gration is possible only if Ar ray < 1, which is violated 
close to the boundary /i = —1. However, in this region 
r ray = r/(l + /x) 3> 1 simply enforces I w S. In the process 
of integration, we set I — S wherever AT ray (/i) > 0.7. 

(2) The transfer equation contains the term (1 — 
a 2 )dl/da = sin# 81/88. We use a grid di (i = 0, ...,n), 
where 9q = and 6 n = n. The term sin 9 81/89 is not 
needed at #o and 9 n (it vanishes). For all other 9i we evaluate 
this term using 81/89 = ili+i ~ h-i)/ ~ 0i-i)- 

The transfer equation (fl2l has no free parameters and the 
solution is unique. The result is shown in Figure 1 . The strik- 
ing feature is the strong beaming of the radiation field in the 
fluid frame, even at large optical depths r ~ 10. Beam- 
ing may be described by the ratio of intensities at fj, = 1 
and /i = —1: b(r) = I(l,r)/I(—l,r). This quantity is 
shown in Figure 2. It significantly deviates from unity start- 
ing at r ~ 10. In the zone of t C 1, 1(1, r) = const and 
J(— 1, r) — Io oc r~ 2 . Therefore b(r) oc r 2 at r ^> R+. 

3.2. Adiabatic cooling 
To examine adiabatic cooling, one should consider the en- 
ergy flux of radiation in the lab frame 7^ = AttIi, where 
h = T 2 [/3(I + I 2 ) + (1 + /3 2 )ii] is the first moment of 
intensity in the lab frameQ Here we cannot take the formal 

4 iirlm are the components of the stress-energy tensor of radiation: 
T 00 = 471-io, T 01 = 4-Kh, and T 11 = 4ttI 2 , where index in T^ v 
corresponds to the time coordinate and index 1 corresponds to the spatial co- 
ordinate in the radial direction. Tensor transformation from the fluid frame 
to the lab frame reads T^ v - 

Aj = Aq = r/3. It gives T 01 



AZA^T^p, where Ag 



A{ = T and 



r 2 [/3(T 00 +T 11 ) + (l + ^ 2 )T 1 
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FIG. 2. — Beaming factor b(r) = 1(1, r)/J(— 1, r) in the fluid frame. 
Dashed curve shows the model with isotropic coherent scattering. Solid curve 
shows the model with Thomson scattering in a cold outflow (Section 4.1). 

limit r — > oo, as the transformation between the lab frame 
and the fluid frame is not well defined in this limit. The total 
luminosity of radiation measured in the lab frame is 

L(r) = 4vrr 2 # = (4tt) 2 r 2 T 2 [0 (I + I 2 ) + (l + p 2 ) h] . 

(14) 

Freely streaming radiation would have L(r) = const. Inter- 
action with the outflow results in adiabatic cooling and L(r) 
decreases with r. 

The adiabatic cooling factor is defined by a(r- ln — > r) = 
L(r) I L(r ln ). We chose r- m <C R* where radiation is nearly 
isotropic; then Ii(r in ) — and ^(nn) = ^o( r in)/3. This 
gives, 



a{r i n ->■ r) = 



L{r in ) 



3 r 2 /3(/q + / 2 ) + (1 + /3 2 )/i 

4 rf n 7 (r in ) 



(15) 

This equation is well-behaved in the limit /3 — >• 1. Adiabatic 
cooling is controlled by Io(r), I\{r), and /2(f), which we 
know from the solution of the transfer equation. Figure 3 
shows the resulting a. In the deep subphotospheric region 
(where radiation is approximately isotropic in the fluid frame) 
a = (V/ri n )~ 2 / 3 as expected. A deviation from this law de- 
velops at t ~ 10 (see also Monte-Carlo simulations in Pe'er 
2008). In the region t <C 1 (r R+) most of radiation 
streams freely and experiences no adiabatic cooling. 

The net effect of adiabatic cooling on the escaping radiation 
is described by 



a(r in -> 00) = 2 [ — 



2/3 



-2rr 2/3 . 



(16) 



3.3. Monte-Carlo simulation 



As an independent check of the results presented above, 
we solved the same transfer problem using the Monte-Carlo 
method. The numerical code is described in B10. In this 
section, we use its simplest version that assumes coherent 
isotropic scattering in the fluid frame. 



FIG. 3. — Adiabatic cooling factor a(rj n — > r). We multiplied a by 
(-R*/ri n ) 2 / 3 in this plot; this combination does not depend on rj n as long as 
r- m <S R*. Dotted line shows the thermodynamic result a = (r/B*)~ 2 / 3 , 
which is valid for approximately isotropic radiation. The actual cooling factor 
obtained from the transfer solution is shown by the solid and dashed curves. 
Dashed curve: model with isotropic coherent scattering. Solid curve: model 
with Thomson scattering in a cold outflow (Section 4.1). 

The Monte-Carlo method is fundamentally different from 
solving the transfer equation. It operates with individual pho- 
tons that are injected at 7'i n <§; i?*. The simulation assumes 
a finite T and is performed in the lab frame. It tracks the 
propagation and random scattering of a large number of in- 
jected photons and accumulates their statistics at different 
radii. These statistics are used to reconstruct the angular dis- 
tribution of radiation intensity in the fluid frame. We ran 
the Monte-Carlo simulation for outflows with T = 20 and 
r = 600. The results were identical, again confirming that 
the transfer does not depend on T as long as T > 1 and 
r(r) = const. 

To reconstruct the radiation intensity at a given radius r 
from the Monte-Carlo simulation, we calculate the angular 
distribution of luminosity passing through the sphere of ra- 
dius r in the lab frame, dL/dft, (r). In the ultra-relativistic 
transfer problem, practically all photons move forward in ra- 
dius and cross a given r only once. We accumulate the statis- 
tics of angles and energies of photons at radius r, which 
gives dL/djl and the intensity of radiation in the lab frame, 
I = (4irr 2 dL / 'dfi. The corresponding intensity in the 
fluid fram e is given by I = T>~ 4 I where T> is the Doppler 
factor (eq. IA3b - Open circles in Figure 1 show the result of 
this calculation at r = R+/2. It is in perfect agreement with 
the solution of the transfer equation. Similar excellent agree- 
ment is found at other radii. 

The strong anisotropy of radiation at subphotospheric radii 
is a result of cooperating effects. Photons with large u have 
a larger free path in the lab frame A (in the opaque zone, 
A ~ [1 + r/r.) The free path is accompanied by a shift in 
H, A/z = (d/.i/ds) A > 0. This leads to the pile up of photons 
along the radial direction. Photons with large u also experi- 
ence less adiabatic cooling than the more frequently scattered 
photons with small u. 
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FIG. 4. — Dist ributio n of the last-scattering radius r*, obtained by inte- 
grating equation )B16t over /j*. Identical distribution is obtained using the 
Monte-Carlo technique. The distribution remains practically the same in 
all models of photon transfer calculated in this paper: isotropic scattering, 
Thomson scattering in cold plasma, and Compton scattering in hot plasma. 
Approximately 2/3 of photons have r* between 0.3/?* and 3R*. 

3.4. Fuzzy photosphere 

R+ was defined in equation (13[ as the radius where the 
parameter r appearing in the transfer equation equals unity. 
It gives an estimate for the characteristic photospheric radius. 
Clearly, the sphere of radius R* cannot be thought of as the 
last-scattering surface, for two reasons: (1) the optical depth 
seen by a photon depends on its emission angle fi = cos 9 
and (2) the free path of a photon near R* is a random variable 
comparable to R±. Therefore, the radius of last scattering is 
a random variable. Its average value logarithmically diverges 
for a steady outflow extending to infinity and cannot be used 
to define the photosphere. 

The process of radiation decoupling from the scattering 
plasma may be described by the probability distribution 
dP / dr±dfi± where /i* = cos^* is the emission angle of the 
photon (measured in the fluid frame) at the last-scattering 
point. Pe'er (2008) considered a similar distribution for r* 
and 6+, where 0* is the emission angle in the lab frame. His 
analytical expression is however inaccurate. The correct ex- 
pression is given in Appendix B. Integrating dP/ dr+d/i+ over 
/z*, one finds the distribution of emitted photons over the ra- 
dius of last scattering, dP/dr± (Fig. 4). 

The same dP/dr* is obtained using Monte-Carlo tech- 
nique. The advantage of the Monte-Carlo simulation is that 
it is easily extended to hot outflows and to scattering with ex- 
act Compton cross section. We calculated dPj d In r* for the 
collisionally heated jet in the model of B10. The result was 
practically identical to that shown in Figure 4. 

In view of the broad distribution of r*, it is not appropri- 
ate to locate the GRB photosphere at any specific radius, as 
emphasized by Pe'er (2008). When a characteristic radius is 
needed for rough estimates, R+ would be a reasonable choice. 
Alternatively, a characteristic photosphere could be defined as 
the sphere outside of which 50% of photons are released (i.e. 
experience the last scattering). The corresponding radius is 



0.8-R*. 

4. TRANSFER WITH EXACT ELECTRON SCATTERING 

Even for cold outflows, the radiative transfer is not ex- 
actly described by the model of coherent isotropic scattering. 
Two effects contribute to this: (1) electron scattering is not 
isotropic, and (2) radiation becomes polarized in the process 
of radiative transfer, with two modes of polarization, so two 
equations describe the transfer problem instead of one. 

Furthermore, if the outflow is strongly heated (as expected 
in GRBs) scattering is not coherent, i.e. does not conserve 
photon energy in the fluid frame. This has a strong effect on 
the intensity of radiation, as instead of passive adiabatic cool- 
ing radiation is heated through the Comptonization process. 

In Section 4. 1 we develop the accurate transfer model for 
Thomson scattering in a cold outflow, which takes into ac- 
count polarization. In Section 4.2 we present a transfer model 
for heated outflows. The results are compared with the simple 
isotropic-scattering model of Section 3. 

4. 1 . Thomson scattering and polarization 

The polarized transfer in a static, cold electron medium was 
described by Chandrasekhar (1960) and Sobolev (1963). In 
axisymmetric problems, e.g. in plane-pallel or spherical ge- 
ometries, there are two polarization modes of radiation: one 
with electric field perpendicular to the plane containing the 
photon direction and the axis of symmetry and the other with 
electric field parallel to this plane. Let I± and Ih be the energy 
intensities of the two modes. Scattering can change the polar- 
ization state, so four scattering processes can occur: _L^_L, 

-L— H|, 1 1 ^ I , and ||— >-||. They are described by four different 

cross-sections. The total intensity of radiation is I = I±_ + I\\ ■ 
The degree of polarization is p = Q/I where Q = I± — In. 

The transfer equations for I and Q in a static medium are 
given in Chandrasekhar (1960) and Sobolev (1963). The gen- 
eralization of these equations to the case of a relativistically 
moving medium is straightforward (see Beloborodov 1998 for 
equations in the plane -parallel geometry). Here we are in- 
terested in spherical ultra-relativistic outflows. We will use 
the corresponding equations for the intensities in fluid frame, 
which are well-behaved in the limit (3 — > 1. The transfer 
equations for I and Q are similar to equation ©El 
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(17) 
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A* 

(18) 

Here S and Q are the source functions in the fluid frame, 
where the scattering medium is static. Their expression in 
terms of moments of / and Q is exactly the same as in the 
static problem described by Chandrasekhar and Sobolev, 

S = I + | (3^ 2 -l)(h-j + Qo- Q2) , (19) 



R = -(1- H 2 



h - j + Qo - Q 2 



(20) 



Note that the polarization states are invariant under Lorentz boosts along 
the axis of symmetry. Therefore, / and Q obey the same transformation 
between the fluid frame and lab frame: I = V 4 I and Q = V 4 Q, where 
T> is the Doppler factor (see Appendix A). The source functions S and R 
transform in the same way. 
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FIG. 5. — (a) Same as Fig. 1 but for Thomson scattering and taking into 
account polarization. The plot shows intensity I(u) (in the fluid frame) at 
three radii: r/R* 0.03, 0.1 and 0.5. (b) The corresponding solution for 
the photon-number intensity X(fi). It is normalized so that X\ + /3Xrj = 1 
(note that X\ [r] + /3Iq [r] = const expresses conservation of photon number, 
see Section 2.3); = 1 in our transfer problem. 

The parameter tt appearing in equations ( TPTT i and (TT~8T > is de- 
fined using Thomson cross section <tt, 



r T (r) = 



O-TiVe 

47rrr 2 (3c 



R* = 



<7TN e 

4nT 2 Pc 



(21) 



where r — R* corresponds to tt = 1. 

Similar to Section 3, we will consider here outflows with 
N e (r) = const and T(r) = const. This implies 7r(r) = 
R*/r and g = 1. Equations ( fTvT > and (fT8T l are solved nu- 
merically in the same way as described in Section 3.1. The 
resulting angular distribution of intensity / in the fluid frame 
is shown in Figure 5a. It is similar to the model with isotropic 



FIG. 6. — Polarization degree of radiation transferred through the outflow 
with Thomson scattering opacity, p = Q/I. Dashed curve shows the polar- 
ization of radiation propagating at angle 8 = tt /2 in the fluid f rame . Solid 
curve shows the angle-averaged polarization defined in equation {22) , 

scattering (Fig. 1). 

The polarization degree of radiation is shown in Fig- 
ure 6. It begins to grow in the subphotospheric region (where 
anisotropy develops), and radiation becomes significantly po- 
larized in the photospheric region. The strongest polariza- 
tion p « 0.6 is found at large radii for radiation at angles 
9 w 7r/2. It is produced by scattering of strongly beamed ra- 
diation, which naturally generates a high polarization at scat- 
tering angles close to 90 degrees. The overall polarization at 
a given radius may be described by the average p that is de- 
fined using the energy fluxes in the two modes measured in 
the lab frame. This definition involves the transformation of 
moments J m and Q m to the lab frame, which gives 



. . Qi 
< p >= — 

h 



fi(Q + Q 2 ) + (l + /3 2 )Qi 
P{h + h) + (1 + P 2 )h 



(22) 



This expression is well-behaved in the limit f3 — > 1 and be- 
comes < p >= (Q + 2Qi + Q 2 )/(Io + 2Ii + h). It reaches 
0.24 outside the photosphere (Fig. 6). 

Equations ( fTTI i and ( fT~8b govern the transfer of energy in 
the two polarization modes. The photon number in the two 
modes is described by intensities Zj_ and Tu . The equations 



fori 
eq.QIB 

dl 



1± + In and Q 



I- 

In read (cf. the similar 



<91nr 

dQ 
<91nr 



(1 



M J 9 



dl 
dfi 

dQ 
Oft 



3 (1 - fig) 2 + t t 



(S-l) 



{I- a 2 ) g-±-3(l-pg)Q + r T 



1 



(23) 

Q) 



i + n 

(24) 

The source functions S and 1Z are related to the moments I m 
and Q m (m — 0, 1, 2) in the same way as S and R are related 
to I m and Q m (eqs.[T9l and 120b. 

We numerically solved equations d23l and d24l i with g = 1. 
The resulting angular distribution I(/x, r) is shown in Fig- 
ure 5b. It is less anisotropic than the energy intensity I (a, r). 
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FIG. 7. — Comparison of photon-number intensity X(fi) at r = R+/2 in 
three models: isotropic coherent scattering (dashed curve), polarized Thom- 
son scattering in a cold outflow (solid curve), and unpolarized Compton scat- 
tering in a hot outflow (open circles), see text for the description of the hot- 
outflow model. In all three models, the intensity is normalized so that the 
conserved quantity X\ + /3Xn equals unity. The models assume /? ~ 1. 

(A similar result is found in the isotropic-scattering model.) 
Correspondingly, the quantity Q /I is smaller than p = Q/I, 
roughly by a factor of ~ 3/4. 

4.2. Scattering in a hot plasma 

GRB jets are heated, which leads to Comptonization — a 
significant flow of heat from particles to radiation. This pro- 
cess depends on the electron temperature that must be cal- 
culated self-consistently. For example, consider the fiducial 
model of the collisionally-heated jet in BIO. The jet with 
Lorentz factor T = 600 is heated at radii r > R n ~ i?*/20, 
and its electron temperature outside R n may be approximated 
by 



= ().()4o I — 



0.23 



(25) 



The initial temperature of radiation at R n is 0.6 keV. Comp- 
tonization begins at R n with Kompaneets' y-parameter y = 
40 e Tr ~ 1- Equation (|23T > is a reasonable approximation in 
the main heating region R+/20 < r < R±. The exact value 
of temperature outside R* is not important — its contribution 
to thermal Comptonization is small, — and the above approx- 
imation for T e will be sufficient for our purposes. 

Comptonization significantly changes the energy intensity 
/(/i, r) compared with the cold-jet model. In addition to the 
thermal plasma, nonthermal particles are continually injected 
with energies ~ 140 MeV from the decay of pions produced 
by nuclear collisions. These particles convert their energy to 
a small number of high-energy photons, which impact the 
radiative transfer and the observed spectrum above 20 MeV 
(B10). 

Here, however, we limit our consideration to the transfer of 
photon number (rather than energy). The high-energy pho- 
tons make a negligible contribution to the number intensity 
I(p, r), so it is sufficient to consider scattering by the thermal 



plasma with temperature (1251 1. which strongly dominates the 
optical depth. The solution for the number intensity l(p, r) 
turns out to be close to that in a cold jet. The obtained angular 
distribution of Comptonized photons at r = R±/2 is shown 
by open circles in Figure 7. Its beaming is somewhat stronger 
compared with the cold-jet model of Section 4. 1, however the 
difference is modest. The effect of electron heating on the 
transfer solution for photon number is ~ 10%. 

5. RADIATION-DOMINATED OUTFLOW 

Anisotropic radiation always tends to push the flow toward 
the equilibrium velocity at which the net flux of radiation 
vanishes in the fluid frame. When this effect is strong, it 
leads to the peculiar "equilibrium transfer" where radiation 
and plasma self-organize to flow with a common velocity (Be- 
loborodov 1998, 1999). In this section, we discuss the condi- 
tions for this regime and the corresponding solution of the 
transfer problem. 

The radiative force applied to each electron (or positron) in 
the fluid frame is given by / = 4ttIi(Tt/c. We assume that 
the electron thermal motion in the fluid frame is slow (non- 
relativistic)0 Then Lorentz transformation gives the same 
value for the force measured in the lab frame, f = f. The 
outflow acceleration is governed by the dynamic equation, 



P 



dim 



0"x 



, = 4ttJi 

ar c 



(26) 



Here I\ is the first moment of radiation in the fluid frame 
(eq. |7]l, p is the rest-mass density of the plasma, and n is the 
number density of e ± ; both p and n are measured in the fluid 
frame. Then, for ultra-relativistic outflows (J3 — > 1), one finds 



<J 



'-St- 



4tt7i 
pc 3 



(27) 



One can use the approximation g = 1 (i.e. T = const) if 
the last term in equation d27l > is much smaller than unity. At 
the photospheric radius R* this term equals AttIi (pc 3 ) -1 = 
x(U/pc 2 ) where \ = h/Io < 1 is a numerical factor that 
depends on the angular distribution of radiation and U — 
47r/ /c is the radiation energy density in the fluid frame. If 
U /pc 2 is not small compared with unity, the transfer equation 
should be solved with the self-consistent function g(r) given 
by equation d27l >. 

The self-consistent solution of the transfer problem can 
be obtained analytically in the extreme radiation-dominated 
regime U/pc 2 1. Note that there exists a special value 



of r = To for which the radiative force (eq. 1261 ) vanishes; 
this value corresponds to Ii = 0. The force is positive if 
r < Tq and negative if T > To, i.e. it always pushes the out- 
flow toward T = To. In the radiation-dominated regime, the 
timescale for dynamical relaxation toward T = Tq is shorter 
than the outflow expansion timescale. This means that the 
outflow maintains T Tq and I\ <C In. ■ 

In this regime, the ultra-relativistic transfer has a simple so- 
lution, 



C 

H^r) = — , 



Q(n,r) = 0, .9 = 0, 



(28) 



6 GRB outflows start with a relativistic temperature kT ~ 1 MeV, but they 
are quickly cooled by adiabatic expansion to much lower temperatures before 
they reach the photospheric radius. Compton scattering provides a strong 
thermal coupling between the plasma and radiation and keeps the electron 
temperature relatively low (non-relativistic). 
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where C is a constant determined by the inner boundary con- 
dition. It is straightforward to verify that equation (f28b is the 
solution of equations ( fTTI i and ( TT8l . The first moment Ix is 
a small (next-order) quantity. It controls the outflow accel- 
eration r a r and is given by I\ = (Attt^)^ 1 pc 3 -c I . 
The outflow accelerates linearly with radius both inside and 
outside the photosphere, as long as U pc 2 . Note that 
U/pc 2 oc r" 1 , so the outflow eventually must reach radii 
where U < pc 2 and the acceleration ends. This transition 
occurs outside R± and has no effect on radiation escaping the 
outflow to distant observers. 

Equation d28l i states that radiation remains isotropic in the 
fluid frame which accelerates as T oc r. The sustained 
isotropy is a consequence of a remarkable fact: a freely prop- 
agating photon between two successive scatterings at radii r\ 
and r2 does not change its angle measured in the fluid frame, 
9i — 02- This fact can be derived as follows. In the lab 
frame the angle of a photon propagating from n to 7-2 satis- 
fies the equation of a straight line, r\ sin Ox = r-i sin 62 ■ Using 
t^/tx = T2/T1 and Doppler transformation sin 9 = Vs'md 
(where V = T[l + cos 9] when /? -> 1), one finds X = 9 2 . 

Thus, free propagation between scatterings does not gen- 
erate any change in the angular distribution of photons 
in the fluid frame, and an initially isotropic radiation re- 
mains isotropic, I = Iq. Scattering of isotropic radiation 
gives isotropic radiation, so the source function also remains 
isotropic in the fluid frame, S = Iq. Naturally, the isotropic 
radiation remains unpolarized. 

It is easy to see why the energy intensity in the fluid frame 
scales with radius as I oc r -4 . First note that conserva- 
tion of photon number implies that the photon-number inten- 
sity scales as I oc r _1 r -2 oc r~ 3 (see Section 2.3 and use 
isotropy, I«2j) and Ix <C Iq). Coherent scattering does not 
affect the photon energy in the fluid frame, so v changes only 
during free propagation of the photon. Propagation between 
successive scatterings at r\ and r 2 occurs with constant en- 
ergy in the lab frame, T>2 hv% = T>x hvx- Using 92 = 9x, one 
finds V2/V1 = Vi/Dx = Tz/T\ = rijrx- Thus, v oc r _1 for 
each photon. Together with I oc r~ 3 this implies I oc r~ 4 . 

The existence of a frame (fluid frame) where radiation re- 
mains isotropic and scattering is coherent implies that scat- 
tering has no effect on radiation, i.e. the transfer occurs as if 
radiation propagat ed i n vacuum. Indeed, consider the basic 
transfer equation (IA1 b in the lab frame. Coherent isotropic 
scattering in the fluid frame gives S u — I v , which implies 
S v = I v in the lab frame and hence dl u /ds = 0. Thus, radia- 
tion intensity in the lab frame remains constant along the ray, 
just like propagation in vacuum. In particular, the spectrum of 
radiation is preserved and its beaming angle decreases as r -1 . 

The corresponding solution for the intensity in the fluid 
frame can be obtained by the Doppler transformation l v = 
T>~ 3 1 V (Appendix A) of the vacuum solution for I v . Alter- 
natively, the same result can be obtained from equation ([TJ. 
Substituting g = 0, one gets 

I v {v,»,r)=I v fc/i,^ (j-^ . (29) 

It confirms that, when viewed in the fluid frame, the trans- 
fer preserves isotropy of photon distribution and shifts each 
photon in frequency as r — 1 . 

To summarize, as long as U ^> pc 2 radiation behaves as 
if there were no scattering and it streamed freely, regardless 



of the optical depth. This is a special feature of the ultra- 
relativistic transfer in spherical geometry. It differs from 
the radiation-dominated transfer with T = To in the plane- 
parallel geometry (Beloborodov 1998; 1999). 

6. VARIABLE JETS AND THE STEADY SPHERICALLY SYMMETRIC 

MODEL 

In this paper, we considered radiative transfer in outflows 
that are steady and spherically symmetric. Here we discuss 
why these assumptions are not so restrictive as they might 
seem and the model may describe variable jets. 

Radial outflows with Lorentz factors r> 1 have two well- 
known features: (1) Their parts are causally disconnected on 
scales larger than l± ~ r/T on any sphere of radius r (r ~ R+ 
should be taken as the characteristic radius for the problem 
of photospheric emission). (2) Since both radiation and fluid 
move outward with almost speed of light, the radial diffusion 
of radiation relative to the fluid is inefficient on scales Sr 3> ill 
where lu ~ r/T 2 . To a first approximation, each "elementary 
pancake" of volume lu x l± X /j_ in the lab frame has its own 
radiative transfer and produces photospheric emission almost 
independently from the neighboring pancakes. If strong in- 
homogeneities of the jet are confined to scales much larger 
than iii and l± (in the radial and transverse directions, respec- 
tively), radiative transfer in each pancake occurs as if it were 
part of a steady, spherically symmetric outflow. 

The independence of emissions from different pancakes can 
be better quantified if one considers the photon exchange be- 
tween two pancakes separated by Sr > In. The exchange is 
one-way only: the trailing pancake can receive photons from 
the leading pancake (the opposite communication is impossi- 
ble for Sr > in). This "trailing diffusion" of radiation was 
studied by Pe'er (2008). He considered a very narrow shell 
of photons (formally a delta-function of radius) injected at a 
small r in a steady jet with a constant Lorentz factor. The pho- 
tons diffuse through the jet and eventually escape, producing 
an isolated pulse of emission that will be received by a dis- 
tant observer. The characteristic observed width of this pulse 
is St bs ~ t* — l\\/c ~ i?*/r 2 c, and it has an extended tail 
whose intensity decreases as (< bsA*)~ 2 ■ This implies that 
the trailing diffusion of radiation on scales Sr > lu is sup- 
pressed as (<5r/Z|| )~ 2 . 

A realistic GRB jet is continually filled with thermal radia- 
tion near the central engine. It may be viewed as a continual 
sequence of elementary pancakes that release their photons 
near R+. The strong Doppler beaming implies that photo- 
spheric emission seen by a distant observer is dominated by 
a small patch l± x l± on the sphere of radius ~ R±, The 
observer receives radiation released by consecutive pancakes 
in the same order as they pass through R+. The observed 
timescale of passage of one elementary pancake through R+ 
is which may be smaller than 1 ms for GRBs. The steady 
transfer model developed in this paper is valid for GRBs with 
variability timescales St Q x, s > t*. The model permits different 
R* for pancakes separated by timescales <5i bs = Sr jc ^S> t*. 

7. DISCUSSION 

This paper explored radiative transfer in ultra-relativistic 
outflows. The transfer problem is well defined and simpli- 
fies in the limit r > 1. In this limit, radiation propagating 
backward in the lab frame can be neglected. Therefore, the 
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transfer solution is independent of the outer boundary condi- 
tion, in contrast to transfer in static media studied by Chan- 
drasekhar (1960) and Sobolev (1963). The problem is solved 
by direct integration of the transfer equation, with no need for 
iterations. The model with T — > oo gives excellent approxi- 
mation to transfer in outflows with finite T > 10. 

7.1. Transfer in radiation- dominated and matter- dominated 

outflows 

The approach developed in this paper gives a simple solu- 
tion for the old problem of radiation-dominated jet discussed 
by Paczyriski (1986) and Goodman (1986). This jet is baryon- 
clean. It is very opaque at small radii because of the thermal 
population of pairs. Almost all pairs annihilate at larger 
radii, and almost all the jet energy is carried by radiation that 
is released at the photosphere. Paczyriski and Goodman con- 
sidered the opaque zone of the radiation-dominated jet and 
derived its Lorentz factor T oc r from energy-momentum con- 
servation. They argued that quasi-thermal emission should be 
observed from the jet, with a spectral peak near 1 MeV. They 
suggested, however, that the observed spectrum should be dif- 
ferent from blackbody because of complicated transfer effects 
near the photosphere. Goodman (1986) performed a numer- 
ical calculation with simplifying assumptions, which gave a 
nonblackbody spectrum. 

In fact, the exact transfer solution for this problem gives 
precisely blackbody spectrum. As shown in Section 5, pho- 
tons in a radiation-dominated jet are transferred as if there 
were no scattering at all. The radiation remains isotropic in 
the fluid frame which accelerates as T oc r both inside and 
outside the photosphere. A distant observer can think that ra- 
diation freely propagates from the central engine of the jetQ 
Note also that regardless of how strong dissipation/heating 
may occur in the jet it does not have a dramatic impact on 
the shape of the spectral peak, because the energy budget of 
heating is negligible compared with the Planck radiation. The 
observed radiation should have a blackbody spectrum. 

The opposite, "matter-dominated" regime was considered 
by Paczyriski (1990). In the opaque zone radiation cools adi- 
abatically and the jet energy becomes dominated by baryons. 
Then it continues inertial expansion (coasting) with some 
relict thermal radiation in it until the radiation is released at 
the photosphere. 

In the matter-dominated regime, the photospheric spectrum 
cannot have the blackbody shape, even if the outflow cools 
passively, with no heating, up to the photosphere. BIO showed 
that the photospheric spectrum in the soft X-ray band has 
the slope a as 0.4 instead of the blackbody (Rayleigh-Jeans) 
slope a = 1. Moreover, collisional heating in GRB jets trans- 
forms the photospheric spectrum into the Band-type radiation, 
with extended high-energy emission instead of the exponen- 
tial cutoff above 1 MeV. Thus, the photospheric spectrum of a 
matter-dominated jet is changed from blackbody both below 
and above the MeV peak. 

Besides giving a non-blackbody spectrum, the radiative 
transfer in matter-dominated jets has other interesting fea- 
tures. Radiation becomes strongly anisotropic in the fluid 

7 The only deviation from the free-propagation solution occurs where the 
jet temperature drops below ~ m e c 2 and the equilibrium density of pairs 
drops below the density of photons. In this region, radiation receives signifi- 
cant energy from the annihilated pairs, which boosts its density by the factor 
of 11/4 (similar to what happens in the expanding universe). After this tran- 
sition the jet is still extremely opaque due to the remaining (exponentially 
reduced) er^ population, and the radiation remains Planckian. 



frame well before it decouples from the fluid. Radiation at the 
characteristic photospheric radius R* has the beaming factor 
b ~ 30 (Fig. 2). Beaming affects the adiabatic cooling of pho- 
tons in the subphotospheric region. The net cooling factor for 
radiation emitted at a radius r- m <C R* equals 2(n n / 'R+) 2 / 3 . 

In a heated jet, adiabatic cooling of radiation is counter- 
balanced (or dominated) by Comptonization, so the mean 
photon energy can grow with radius. This has a strong effect 
on the transfer solution for the radiation intensity. However, 
if one focuses on the transfer of photon number (rather than 
energy), the results are not sensitive to heating. The photon- 
number intensity I in the passively cooling and heated jets is 
very similar (the difference is ~ 10%). In both cases, X is 
strongly beamed in the subphotospheric region (Fig. 7). 

7.2. Detecting the blackbody component in GRBs 

As discussed in Section 7.1, photospheric emission in 
GRBs can have a blackbody spectrum only when the photo- 
sphere is dominated by radiation (i.e. U ^> pc 2 at r ~ i?*). 
The detection of a blackbody component in a GRB spectrum 
would provide clear evidence that part of the photospheric 
emission is in the radiation-dominated regime. 

The existing data are inconclusive. It includes 
GRB 090902B that was much discussed recently as a 
burst with a blackbody component. In fact, it is equally well 
fitted by the Band function plus a power law (Ryde et al. 

2010) . The data interpretation is further complicated by the 
variability of photospheric emission. It may vary on very 
short timescales, as short as t+ ~ R+/cT 2 , which can be 
smaller than one millisecond for a typical GRB. The achieved 
temporal resolution of spectral analysis is far worse than 
1 ms and may not give the true instantaneous photospheric 
spectrum. The photospheric emission may quickly switch 
between the radiation-dominated and matter-dominated 
regimes and these variations would remain undetected. 

The low-energy slopes of the observed GRB spectra, a, 
are affected by the time averaging, which tends to reduce a. 
Bursts with largest a are most promising for detecting the 
blackbody component. In some cases, a ~ 1 were reported 
(Ghirlanda, Celotti & Ghisellini 2003). This indicates the ex- 
istence of the radiation-dominated regime. 

7.3. Detecting photospheric polarization 

Radiation remains unpolarized in radiation-dominated jets 
(Section 5). In contrast, in matter-dominated jets, radiation 
acquires a strong linear polarization in the photospheric re- 
gion (Section 4.1). An ideal detector that has enough angu- 
lar resolution to image the spherically-symmetric outflow on 
scales ~ r _1 i?* would detect the polarization. In practice, 
such a high angular resolution is not achieved. The detectors 
receive a mixture of radiation whose polarization averages to 
zero unless something breaks spherical symmetry. 

Three principle possibilities for breaking the symmetry are 
as follows: (1) The main emitting region of size l± ~ r _1 i?* 
is partially eclipsed. (2) The outflow deviates from spherical 
symmetry on scales ~ l ± , (3) The jet carries magnetic fields 
with a coherence scale l±. In magnetized jets, the syn- 
chrotron component of photospheric emission becomes dom- 
inant at photon energies below ~ 100 keV (B10; Vurm et al. 

2011) . This component can be highly polarized. Future po- 
larization measurements across the X-ray spectrum will help 
estimate the magnetization of GRB jets. 
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7.4. Modeling frequency-dependent radiative transfer in 
heated jets 

The photosphere of a baryonic jet is a fuzzy object — about 
2/3 of photons are released in the region i?*/3 < r < 3i?*, 
and the remaining 1/3 comes from even more extended region. 
Modeling the heated anisotropic radiation emerging from this 
region requires accurate transfer simulations. 

BIO developed a Monte-Carlo transfer code that solves the 
transfer problem in a broad range of photon energies up to 
100 GeV, including the effects of 7-7 absorption. Alterna- 
tively, one can use the kinetic method that solves the kinetic 
equations for the photon and electron distribution functions 



(Pe'er & Waxman 2005; Vurm & Poutanen 2009). The devel- 
oped kinetic codes have, however, one drawback: they assume 
isotropic radiation in the fluid frame, which is not a good ap- 
proximation. Besides, it violates conservation of photon num- 
ber in the lab frame. The kinetic method can be used more 
efficiently if it calculates the evolution of radiation by solving 
the transfer equation (Q~|i. This method will be implemented in 
an upcoming paper (Vurm et al. 201 1). 



This work was supported by NSF grant AST-1008334 and 
NASA grant NNX10AO58G. 



APPENDIX 

A. BASIC EQUATIONS OF RELATTVISTTC TRANSFER 

Consider radiation with specific intensity /„ (i>, /i, r) in the fixed lab frame. Here v is the photon frequency, fl — cos 8, and 8 is 
the photon angle with respect to the radial direction. Hereafter quantities measured in the fixed lab frame are denoted with tilde. 
The transfer equation reads 



as 



I), 



(Al) 



where ds is the path element along the ray and h v (v, fi, r) is the absorption coefficient of the outflow in the lab frame. The source 
function S v equals j v /K v , the ratio of emission and absorption coefficients (e.g. Chandrasekhar 1960). 

The outflow in our problem is moving radially with velocity f3(r) and Lorentz factor T(r). The transfer equation in the lab 
frame is not well behaved in the ultra-relativistic limit r > 1. Therefore, we rewrite it in terms of intensity J„ measured in 
the fluid frame, i.e. in the frame comoving with the outflow. This is straightforward to do using the usual transformation laws 
(Prokof'ev 1962; Castor 1972; Mihalas 1980). The transformations are given by 



2V, 



= 23 3 /„ S v = 
V = T(l-j3fx)- 
- ^ + /? 



V 3 S„, 

1 =r(i + /3/x), 



Using 7„/z> 3 = I v /v 3 and S v /v 



r,3 



1 + /3/X 

S v jv z \ equation dAU may be written as 



l v is considered as a function of r, /1, v, and the derivative along the ray is expanded as 

dl v dr dl v d[i dl v dlnv dl v 
ds dr 



ds 



ds dji 



ds d In v 



(A2) 
(A3) 

(A4) 
(A5) 

(A6) 
(A7) 



Here one can use dr/ds = jl, dv/ds = 0, and djl/ds = (1 — jl 2 )/r (a consequence of r sin 8 = c onst , wh ich is valid for any 
straight line). The corresponding derivatives of v and /1 are obtained using the transformations (IA3b and ( IA4b . This gives 



dr 
ds 
dfj, 
ds 
din v 
ds 



r 

v 

1 - 



2 



1 



r . a . dinr 
vp {fi + p) dh^ 



TP (l-/x 2 ) 
V r 



I>. Q . rflnT 



where we used dj3 = dr//3T 3 and (1 — p 2 ) =T> 2 (1 — /i 2 ). Then equation JA6t becomes, 



(J*+P) 



where 



dlv 
dlnr 



1 + 



(/i + /3) dlnT 



dlnr 



^(1 



P 



dhir 



d In v 



3/. 



(A8) 
(A9) 

(A10) 



(All) 
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Equation (2.12) in Mihalas (1980) is reduced to equation (IA1 lb in the steady case. 

B. ANALYTIC SOLUTION 

In Section 2 we solved the transfer equation numerically. Here we collect useful analytical formulas that may be used instead 
of the numerical solution. 

Optical depth along the ray 

Consider a photon propagating from radius n to T2 along a straight line in the lab frame. Let di be photon angle at n . The 
optical depth along the ray from r\ to r-i is 



Tray^lA,^) 



re(r,0) 



dr 
cos 9 



(Bl) 



Here 6{r) is photon angle at radius r. It satisfies the relation (which expresses the fact the photon moves along a straight line), 

r sin Q = r\ sin 9\ . (B2) 

The scattering opacity in the lab frame is given by k = T>~ 1 k = T(l — j3 cos 9) an. Let us consider an outflow with T(r) = const 
and n oc r~ 2 . Then the elementary integral in equation ( IB II ) gives 



Tray(ri,0i,r 2 ) = T(ri)T 2 



sin 9\ 



■0(1-^ 

f2 



where r(r) = nar/T. If T 1, one can expand equation (IB3b in T 1 , 



rrayCrijXijra) = 



6 



1-H 

7*2 



3 



1 + — + 4 «i 



'2 



O T- 



Here a; = 0(1) is a convenient variable related to the photon angle 

x = r 2 e 2 = 

1 + /x 

where /x = cos 9 is measured in the fluid frame. 
If we take T2 — > oo and drop index "1" for the emission point, equations (IB3b and (1B41 > are reduced to 



r ray (r, x, oo) = r(r) T 2 



sin 6> 



/3 



-(r) 
6 



(3 + x) + e»(r- 



(B3) 



(B4) 



(B5) 



(B6) 



A similar formula for the optical depth along the ray from radius r to infinity is given in Abramowicz et al. (1991) and Pe'er 
(2008). 

Expressions for intensity and source function 

The formal solution for the transfer problem is written in terms of the source function. Let us first consider the transfer of 
photon number (Section 2.3). The corresponding intensity in the lab frame is given by 



T(p,,r) 



5(/ii,n) exp[-T ray (ri,/ii,r)] eJT ray (ri,/2i,r), 



(B7) 



where r\ < r is running along the ray and /xi = cos 9\ is related to fi = cos 9 by equation dB2b . Equation jBlj can be rewritten 
in terms of r) and <S(/i, r) using the transformations 2 = T> 3 I and 5 = V 3 S. For outflows with /3 — > 1 this gives 

3 



X(fi, r) 



1 + AH 
1 + A* 



exp [-r ray (n , Mi, r)] dr ray (n, ^i , r). 



(B8) 



Here one can substitute equation ( IB4l i for r ray (ri, /xi, r). The identity r ray (ri, /xi, 



-T ray (r, /x, ri) simplifies the integral. 



It is sufficient to know the source function 5(/x, r) to reconstruct the solution for X(/x, r). In the model with isotropic scattering 
the source function S — T$ (r) does not depend on zx. Our numerical result for Xq (r) agrees within a few percent with the 
following formula, 



r A Xq (r) = JC I — H — arctan 

2 7T 



l /-R* 



(B9) 



where /C is a constant. It can be expressed in terms of the photon flux in the lab frame T (Section 2.3), which satisfies r 2 T[r) 
const. At r/R± > 1 we have X\ = X and T = 4t7T(1 + /3)X . This gives (with /3 -> 1) 

r 2 J-(r) 



8tiT 



(BIO) 
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Substitution of equation JB9b to equation 
integration of the transfer equation. 
Similarly, the formal solution for transfer of energy is given by 



and its integration over radius recovers I(/x, r) that was found by direct numerical 



1 



Mi 



exp | 



'ray 



The model with isotropic coherent scattering has S(u, r) = Io(r). In this case, one can use the following formula, 



r 2 I (r)=K [- + 



2/3 



(BID 



(B12) 



where constant K is determined by the inner boundary condition at r <C i?*. Equation (IB 12b remains approximately valid 
for cold outflows with Thomson scattering. For heated jets with significant Comptonization of radiation, Io(r) and S(u, r) are 
different and depend o n the heating history. 

In contrast, equation JB9b remains an excellent approximation even for heated outflows, as long as scattering is the main source 
of opacity and emissivity. 

Distribution of the last-scattering radius and angle 

Consider all photons escaping to infinity from a steady, spherically symmetric outflow. Let N be the number of escaping 
photons per unit time, measured in the lab frame. One can think of N as the rate of photon emission by a source distributed 
throughout the volume of the outflow and attenuated by the optical depth. The photon emission rate from volume element dV 
into solid angle dCl is e dV dVt, where e(r*, 6*) is the photon emissivity in the lab frame. It depends on the radial position of the 

emitter dV, r*, and the emission angle with respect to the radial direction, 9+. The attenuation factor is exp[— T ra y(r*, 8+, oo)], 
which gives 

rray(f*,0*,OO) 



dVdQ 



e exp 



(B13) 



A similar expression holds for the angular distribution measured in the fluid frame (note that e d(l = e d£l is invariant under 
Lorentz transformation). Substituting dV = Airr 2 dr+, d£l = d(f)± d(J,+ (where //* = cos 6>* describes the emission angle in the 
fluid frame), and integrating over 0*, we obtain 



dN 



= 8n 2 rl e exp [— r ray (r*, /i*, oo)] 



(B14) 



dr+du+ 

This equation describes the distribution of escaping photons over the last-scattering radius and angle. The distribution can be 
normalized to unity if we divide it by N = ^r^Tir*) where T is the photon flux in the lab frame. This gives the probability 
distribution for r* and 

dP S r . 

- — - — = 27TK— exp -T ra y(r*,/x*, oo) . (B15) 
dr* d^ J 

Here we used the relation e = kS in the fluid frame, where n — an is the scattering opacity. Equation ( IB 15b shows that the 
distribution of the last-scattering radius and angle is proportional to the source function in the fluid frame, S, whi ch is dete rmin ed 
by the tra nsfer solution. For the isotropic-scattering model one should use S — Io(r). Substitution of equations ( IB6I ) and ( lB9b to 
equation ( IB 15b gives 



dP 



d In T* dfj,+ 4r* 



• arctan 



R~* 



exp 



6R* 



(B16) 



The corresponding distribution of r* and p,* is given by dP / dr+dfl* — T> 2 dP/dr+d/i*. 
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